
#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <time.h>

/*
__global__ void SimpleCorrValue(unsigned char* actual, unsigned char *template, 
        int width, int height, int n_live_pixels, float* result, size_t pitchActual, size_t pitchTemplate) {
/**
 * ASSUMPTION: data is of even length 
 */
__global__ void mean(unsigned char* data, int length, unsigned char* mean) 
{
    
    // Load block area in data to the shared memory
    const int blockSize = blockDim.x;
    __shared__ unsigned char blockData[blockSize];
    int threadPos = blockDim.x * blockIdx.x + threadIdx.x;
    blockData[threadIdx.x] = data[threadPos];
    __syncthreads();
    // Fold the shared array
    int numThreads = blockDim.x;
    while(numThreads > 1) {
         numThreads = numThreads / 2;
        blockData[threadIdx.x] += blockData[threadIdx.x + numThreads];
       
        __syncthreads();
    }
    
    
    // Mean is now on blockData[0]
    // Move mean to the global heap
    const int gridSize = gridDim.x;
    __device__ unsigned char globalData[gridSize];
    globalData[blockIdx.x] = blockData[0];
    __syncthreads();
    // First thread in each block folds the global heap.
    if(blockIdx.x == 0) 
        numThreads = gridDim.x / 2;
        while(numThreads > threadPos && numThreads > 1) {
            globalData[blockIdx.x] += globalData[blockIdx.x + numThreads];
            numThreads = numThreads / 2;
                __syncthreads();
        }
    __syncthreads();
    
    // Divide by length and put globalData[0] into "mean".
    if(blockIdx.x == 0) *mean = (globalData[0] / length);
}
/*
__global__ void standardDiv(unsigned char* data, int length, unsigned char* mean, 
        unsigned char* standardDiv)
{
    // Load block area in data to the shared memory
    __shared__ unsigned char blockData[blockDim];
    int threadPos = blockDim.x * blockIdx.x + threadIdx.x;
    blockData[threadIdx.x] = data[threadPos];
    __syncthreads();
    // Subtract the mean from every data member
    blockdata[threadIdx.x] -= mean;
    // Square every data member
    blockdata[threadIdx.x] *= blockdata[threadIdx.x];
    __syncthreads();
    // Fold
        int numThreads = threadDim.x / 2;
    while(numThreads > threadPos && numThreads > 1) {
        blockData[threadId.x] += blockData[threadId.x + numThreads];
        numThreads = numThreads / 2;
        __syncthreads();
    }
    // Now the sum is at blockData[0]
    // Sum goes to the global heap
    __global__ unsigned char globalData[gridDim];
    globalData[blockIdx.x] = blockData[0];
     __syncthreads();
    // First thread of every block folds the global heap
     if(blockIdx.x == 0) 
        numThreads = gridDim.x / 2;
        while(numThreads > threadPos && numThreads > 1) {
            globalData[blockIdx.x] += globalData[blockIdx.x + numThreads];
            numThreads = numThreads / 2;
                __syncthreads();
        }
    __syncthreads();
    // divide by length and put globalData[0] into standardDiv
    if(blockIdx.x == 0) standardDiv = (globalData[0] / length);
}
*/
int main(int argc, char* argv[]) {
    unsigned char* sampleData = (unsigned char*)malloc(sizeof(unsigned char) * 6);
    
    sampleData[0] = 55;
    sampleData[1] = 55;
    sampleData[2] = 55;
    sampleData[3] = 55;
    sampleData[4] = 55;
    sampleData[5] = 55;

    
    unsigned char* deviceData;
    cudaMalloc(&deviceData, sizeof(unsigned char) * 6);
    cudaMemcpy(deviceData, sampleData, sizeof(unsigned char) * 5, cudaMemcpyHostToDevice);
    unsigned char* devMean;
    cudaMalloc(&devMean, sizeof(unsigned char));
    const int blocksNum = 2;
    const int threadsNum = 4;
    const int length = 6;
    mean<<<blocksNum, threadsNum>>>(deviceData, length, devMean);
    unsigned char* mean = (unsigned char*)malloc(sizeof(unsigned char) * 6);
    cudaMemcpy(devMean, mean, sizeof(unsigned char), cudaMemcpyDeviceToHost);
    printf("Mean is %d\n", mean);
    return 0;
}